Non-Abelian statistics as a Berry phase in exactly 
solvable models 



Ville Lahtinen and Jiannis K Pachos 

School of Physics & Astronomy, University of Leeds, Leeds LS2 9JT, UK 
E-mail: pyvtl@leeds.ac.uk (Ville Lahtinen) 

Abstract. 

We demonstrate how to directly study non-Abelian statistics for a wide class 
of exactly solvable many-body quantum systems. By employing exact eigenstates 
to simulate the adiabatic transport of a model's quasiparticles, the resulting Berry 
phase provides a direct demonstration of their non-Abelian statistics. We apply this 
technique to Kitaev's honeycomb lattice model and explicitly demonstrate the existence 
of non-Abelian Ising anyons confirming the previous conjectures. Finally, we present 
the manipulations needed to transport and detect the statistics of these quasiparticles 
in the laboratory. Various physically realistic system sizes are considered and exact 
predictions for such experiments are provided. 
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1. Introduction 

A striking feature of topological phases of matter is that they can support anyons. 
These are quasiparticles with statistics different from bosons or fermions. The statistical 
behavior of anyons is demonstrated by their adiabatic exchange, which causes non-trivial 
evolution in their quantum state. For Abelian anyons the evolution is given by a phase 
factor. The presence of non-Abelian anyons gives rise to degeneracy in the energy 
spectrum and the evolution is described by a unitary matrix acting on the degenerate 
states. In general, anyons are known to exist in different varieties distinguished by 
their characteristic statistics pQ. Therefore, the explicit demonstration of the statistical 
behavior is essential for the unique characterization of a topological phase. Such phases 
are of great interest due to the possibility of realizing anyons in a physical system 
and due to their potential for technological applications. In particular, topological 
quantum computation employs the anyonic statistics for performing error-free quantum 
information processing [2]. 

The best known many-body system conjectured to support non-Abelian statistics 
is the fractional quantum Hall liquid 0, HI IS]. Other proposals include the p-wave 
superconductor j5j E] as well as various lattice models [3, El E]. These systems are 
either tailored to identically support non-Abelian statistics and have complex physical 
realizations, or they can be described by simple Hamiltonians, but their statistical 
behavior is based on indirect arguments. In particular, for the fractional quantum 
Hall states they rely on properties of trial wave functions [TO], [TTj, [12], whereas for the 
lattice models explicit calculations have not been previously attempted. Although the 
indirect arguments are sound, direct calculations of the statistics are crucial to resolve 
any ambiguities, to address physical realizable finite-size systems and to provide exact 
predictions for the experiments. 

Here we demonstrate how to directly calculate the non-Abelian statistics for a 
class of exactly solvable models. By applying the Berry phase technique [10] to the 
Kitaev's honeycomb spin lattice model [9], we calculate the evolution associated with 
an adiabatic exchange of quasiparticles. This is performed using exact eigenstates of 
a 360 spin system. We obtain a unitary matrix that corresponds to the statistics of 
the conjectured non-Abelian Ising anyons. Together with the fusion rules of these 
anyons [131 EH]) this conclusively demonstrates the non-Abelian character of Kitaev's 
model, thereby confirming the conjectured behavior. Further, we present a scheme for 
creating, transporting and characterizing the anyons that could be used in the proposed 
physical implementations [15J and provide exact predictions for a physically realistic 
range of the model's parameters. 
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Figure 1. (a) The honeycomb lattice on a torus containing two vortex pairs. This 
vortex configuration is created by setting Uij = -1 on the links crossed by solid lines 
and — 1 on all other links. The parameter d controls the minimal vortex separation. 
It is related to the torus dimensions through M = 2 (2d + 4) and N = 2d + 3 (picture 
not on scale). The four dashed arrows Ci, Cf 1 , C 2 and C^ 1 are the oriented parts of 
the path C along which the vortices are moved, (b) C; = C\C%Cy C^ 1 is topologically 
equivalent to a link, (c) C — CiC-f 1 C^C.^ 1 is topologically equivalent to two unlinked 
loops. 

2. The honeycomb lattice model 

Kitaev's model [9] comprises of spin-1/2 particles residing on the vertices of a honeycomb 
lattice. The spins interact according to the Hamiltonian 

H = ~ E ~ E KmrfOjOt, (!) 

u£{x,y,z} £v-links 

where are positive nearest neighbor couplings on links of type v (see Figure[]Ja) 
for link labeling). The second term is an effective magnetic field with positive next-to- 
nearest neighbor couplings such that every plaquette p contributes the six terms 

(i,j,k)ep 

The enumeration of the sites is shown in Figure Ufa). The Hamiltonian has the 
symmetry [H,w p ] = 0, where w p = afa^c^afa^crQ are plaquette operators whose 
eigenvalues w p — — 1 are interpreted as having a vortex on plaquette p. We represent 
the spin operators as a\ = ib^Ci, where Cj, bf,b^ and b\ are Majorana fermions [91 [13]. 
Subsequently, the Hamiltonian takes the form H — | £\ . A^CjCj, where 

Aij = 2Jijiiij + 2 K ijk u ik u jk , 11^ = ibib u j. (2) 

Here Jjj and it^- are shorthand notations for J^- and u\j when is an z/-link. Since 
the mapping to Majorana fermions doubles the size of the Hilbert space, the eigenstates 
of the original Hamiltonian ([I]) are subject to the constraint 

A I *) = !*>, D = [A,<]=0, (3) 
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which follows from the operator identity cfa\af = bfb^b^Ci = 1. Since [H,v,ij] = 0, the 
Hilbert space splits into sectors each labelled by u, a certain pattern of eigenvalues Uij. 
The configurations u can be understood as a classical Z 2 gauge field with local gauge 
transformation operators Dj. Consequently, the plaquette operators w p = Y[(ij)ep^v 
can be identified with gauge invariant Wilson loop operators, whose patterns of 
eigenvalues label the physical sectors of the model. Fixing the gauge field configuration 
u gives then a particular vortex configuration. Throughout this paper we use J u , K and 
u without indices to denote global configurations of these local quantities and use indices 
only when referring to their local values. For instance, K = a means that = a for 
all i, j and k. 

Consider the system of 2MN spins on a torus and assume u to be fixed such that 
it creates the four vortex configuration shown in Figure [T](a). Diagonalization reduces 
the Hamiltonian to the canonical form H = Y2k=i e k\b\b k — ~], where bk are fermionic 
operators satisfying {bl,b t } = 5^ and are the corresponding positive eigenvalues. In 
[T3] it was shown that when the system is in the non-Abelian phase (J x — J y — J z — 1, 
K > 0), the presence of 2n well separated vortices gives rise to n zero modes (e^ ~ 0, for 
k = 1, . . . , n). Importantly, these are separated from the rest of the fermionic spectrum 
by a finite energy gap. In our case of four vortices this implies fourfold ground state 
degeneracy arising from a pair of zero modes that can be either occupied or empty 

i*a 1Q2 > =(&ir(4rigs>, (4) 

where ai,a2 = 0,1 and | gs) = Y\!k=\ ^ I 4) * s ^ ne ground state. For convenience we 
choose the reference state such that b* k I <t>) — 0. 



Numerical diagonalization of A, (see ([2])), gives 2MN eigenvectors ip k satisfying 
the double spectrum Aip^ = iekip^, where coincide with the positive eigenvalues 
of the diagonalized Hamiltonian. We construct a representation of the two degenerate 
ground states | ^10) and | ^oi) as 

MN-l 
k,...,l^a} 

where a = 1, 2, respectively, and Sk,...,i is the fully anti-symmetric tensor of rank MN — 1. 
In general, such states are too large to be stored in a computer, because their number of 
elements grows exponentially with the system size. However, the inner product of two 
such vectors, each depending possibly on some parameters t and t', can be efficiently 
calculated and is given by 

(* a (t)|* /J (0>=det(s5), (6) 

where [B%] kl = V* f (*)Vf (*')• 

2.1. The Ising any on model 

It has been conjectured that Kitaev's model supports the Ising anyon model [91 [3j [5]. 
This model has three types of particles: 1 (vacuum), ip (fermion) and a (non-Abelian 
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anyon). In [12] these are identified with the ground state, the fermion modes and the 
vortices, respectively. The non-trivial fusion rules are given by ip x ip — 1, ip x a = a 
and a x a = 1 + tp. The last fusion rule implies that there is a degree of freedom 
associated with the different ways a number of cr's can fuse when their total anyonic 
charge is fixed. Taking the four a particles to fuse to a ip, this fusion degree of freedom 
is encoded in the two dimensional fusion space, V^ 4 . Its basis can be chosen to be the 
states associated with the two distinct pair-wise fusion channels: 

(cr X cr) x (a x cr) — > tp x 1 = ■0, . . 

(cr x a) x (a x a) — > 1 x tp = ip. 

In [13] the number of intermediate ip's is identified with the number of occupied zero 
modes. Hence, a suitable basis is given by the states {| , | ^2)} (see ([5])). The braid 
operator, R, describes the statistics of the cr anyons. In particular, the monodromy 
operator, R 2 , corresponds to one particle encircling another clockwise. On the basis (GO) 
the monodromy of two u's that belong to different pairs is given by 




3. Non-Abelian statistics as a holonomy 

When z\ and Z2 are the coordinates of the a anyons, their statistics is given by the 
transformation of the wave function under their permutation, i.e. ip(zi, z 2 ) = Uip(z2, zi) 
with U being the characteristic statistical phase or matrix. In real physical systems the 
permutation of the coordinates corresponds to adiabatically transporting the anyons 
such that their positions are swapped. When the positions are swapped twice, i.e. a 
particle winds around the other along a suitable chosen path, the statistics corresponds 
to the accumulated wave function evolution, which is given by the Berry phase, or 
the holonomy [T0l[T6]. For bosons and fermions this is always trivial, with non-trivial 
evolution being a sign of anyonic statistics. 

We demonstrate the statistics of a anyons by adiabatically transporting a vortex 
around another. Consider a Hamiltonian H(X) with n-fold degeneracy {| ty a (\)) \a = 
l,...,n} that depends on some parameters A. When we adiabatically vary A along 
a closed path C, the evolution of the degenerate subspace is given by the holonomy 
r c = Pe^§ c A^(\)d\^ where [A^(\)] a(3 = (* a (A) | I ^(A)) and P denotes path 
ordering in A. To simulate the vortex transport, we discretize the path C into T 
infinitesimal intervals of length SX with A(t) denoting the control parameter value at 
step t. It follows that the holonomy takes the form 

r ^= T i ^ p n (i>«(^))><^( A w)i) , (9) 

t=l \a=l / 

i.e. in the limit 5X — * it is given by the ordered product of projectors onto the ground 
state space at each step t. 
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We evaluate the evolution in the fusion space V%. The basis states <$5§ are not 
symmetrized under gauge transformations ([3]). Nevertheless, their holonomy coincides 
with the holononomy of symmetrized states when C is a loop in both the space of four 
vortex and gauge field configurations. This is due to the orthogonality of states belonging 
to different sectors of u. A suitable path is illustrated in Figure [0(a), where the path C 
(dashed lines) is split into four parts. Different ordering of these parts corresponds to 
the topologically inequivalent paths C\ and C Q given in Figure [Jib) and (c), respectively. 
Neither path spans any area and hence all contribution to the holonomy is topological. 
Since u is a static background field, we need to introduce classical control parameters to 
physically implement the transport. Assuming local control of and on all links, 
we see from (T5]) that the simultaneous sign change of these quantities on link (i,j) is 
equivalent to changing uij — > — Uij. This either generates a vortex pair or transports a 
vortex through the link (i, j). In our simulation this is performed in S infinitesimal steps. 
Taking A = (J, K) and assuming T to be a sufficiently large, the discrete holonomy Q 
for the degenerate states (jSJ) is well approximated by 

p jj( det(^ +1 ) det(i^ +1 ) \ 
^11 1 det(i4 m ) det(£^ I ' 



c 

t=i 



22 



where we have used the inner product (EjLjj. Therefore, the holonomy can be evaluated 
by diagonalizing the Hamiltonian at each step t and multiplying together the inner 
products of the eigenstates from successive steps according to ( |T0|) . We perform this for 
the three parametrizations shown in Table HJ 



Table 1. Three parametrizations (i), (ii) and (iii) for which the holonomy is evaluated. 
Here T = 8S{d + 1) and the number of spins is 2MN = 8(d + 2) (2d + 3). S has been 
increased in (iii) to suppress accumulation of discretization errors due to longer path. 





d 


S 


T 


2MN 


(i) 


1 


2-10 3 


32-10 3 


120 


(ii) 


2 


2-10 3 


48-10 3 


224 


(iii) 


3 


4-10 3 


128-10 3 


360 



Since the spectrum varies slightly with t during the braiding process, we define the 
minimal fermion gap, A, and the maximum energy splitting between the two ground 

states, 5, by 

A = min(e* - 4), 6 = max(e* - e\), (11) 

respectively, where e\ is the fcth eigenvalue at step t. These are plotted in Figure El 
where we observe that both the fermion gap and the level of degeneracy improve as 
K and d increase. Under the adiabatic approximation the holonomy corresponds to 
the exact time evolution when A ^> 5 and 5 0. To physically accommodate these 

| The freedom in changing the basis at each step t of the path C gives rise to an accumulated unitary 
matrix, M, making the result of the adiabatic evolution to be in general given by MTc [12]. Here we 
choose | ^ a (t = 0)) = | V a (t = T)) which gives M = 1. 
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Figure 2. The minimal fcrmion gap A ( ) and the maximum energy splitting 

between the ground states S, (fTTj) (----) as functions of K for parametrizations (i) 
(O ), (ii) (□) and (iii) (()) given in Table Q] The fermion gap grows linearly and the 
degeneracy improves with increasing K for all parametrizations. The fermion gap is 
relatively insensitive to the vortex separation, whereas the degeneracy improves when 
the vortices are further apart. 



conditions in a finite size system, the vortex transport should be fast enough compared 
to 5 for the states | ^ a ), a = 1, 2 to appear as degenerate, but slow enough compared 
to A so that no fermionic excitation is produced. We see from Figure [2] that 4 < 10~ 2 
for (iii) when K > 0.07. This region can support the adiabaticity conditions and hence 
we take K « 0.07 as a lower bound for identifying a stable topological phase. 

To quantitatively study the holonomy, we introduce a fidelity measure for a target 
matrix U and a test matrix V as 

s(U,V) = -tr(UV* + Vlfl). (12) 

For U, V unitary 2x2 matrices we have that s(U, V) — 1 if and only if U — V, 
while in general s(U,V) < 1. Let be the numerically obtained holonomy with 
off-diagonal elements re* e , < r < 1. After fixing the gaug^U, we evaluate the 
unitarity measure, s(l, r^r^ ), Figure [3]^a), and the two different fidelity measures 
of the holonomy: s(|.R 2 |, (TcJ) = r (measure of off-diagonality that characterizes R 2 ) 
and s(R 2 ,T Cl ) = \[s{R 2 , T Cl ) + 1] = |[rcos(f + 0) + l] (the total fidelity), Figure E£b-d). 
Here \U\ denotes a matrix U with its elements replaced by their absolute values. 

First, we observe that the unitarity measure is above 98% for all parametrizations 
when K < 0.10, which we take as an upper bound for identifying a stable topological 
phase. For (i) we obtain no significant off-diagonality due to the small size of the 
system. However, for (ii) the holonomy is predominantly off-diagonal (e.g. r > 0.9) for 
0.02 < K < 0.04, and for (iii) for 0.02 < K < 0.09. The total fidelity, s, accounts also 
for the overall phase and can distinguish between the Ising (s = 1) and SU{2) 2 {s = \) 

§ The holonomy (JTUJ) is only given up to a gauge transformation g : Tc — > g^c9 [IS]- Before s(U, V) 
can be evaluated, the gauge g must be fixed. Due to the finite size of the system the two ground states 
are never perfectly degenerate (see Figure [2]), implying g = diag(e Ic ^ >1 , e 1 ^ 2 ) for some random phases 4>i 
and 4>2- This can be easily taken into account. 
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Figure 3. (a) The unitarity measure, s(l, r^r^ ), as a function of K for the three 

configurations given in Table [TJ The measure of off-diagonality, s(|i2 2 |, \^c t \) ( )> 

and the total fidelity, s(i? 2 , Tc,) (- - - -), as a function of K for the parametrizations 
(b) (i) (O ) ; (c) (ii) (□) and (d) (iii) (0). Based on unitarity and the energy gap 
behavior, we expect a stable phase in the area 0.07 < K < 0.10 bounded by the 
dashed vertical lines. 

anyon models whose monodromies only differ by an overall phase factor, e~ ln ^ 2 pp. We 
observe that for (ii) there is a small region around K 0.02 and for (iii) there is a wider 
region, 0.08 < K < 0.10, where s > 0.9. The maximum fidelities are given by 0.981 
and 0.991, respectively. Parametrization (iii) also has a region 0.02 < K < 0.05 where 
s ~ | with error ±10 _1 . However, we disregard this regime, because such a region does 
not exist for the smaller system (ii) and it lies outside the domain which we consider 
as a stable topological phase. Further, we check for all parametrizations and all K that 
Tc ~ 1 with error less than 10~ 2 , that for K — the holonomy vanishes and that 
T c -i = when the direction of braiding is reversed. 
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4. Braiding and detection of the non-Abelian statistics in a laboratory 

Since our calculation involves only the experimentally accessible parameters J and K, 
it translates directly to how one could physically implement the creation and transport 
of anyons in the laboratory. In particular, in the optical lattice proposal of Micheli 
et al. [E], the vortex transport would correspond, given sufficient site addressability, 
to the local adiabatic inversion of magnetic field as well as of the couplings J through 
the introduction of suitably tuned lasers. Also, as the energy of the system is known to 
depend on the zero mode populations (13] , the effect of braiding can be detected through 
spectral means. As the monodromy swaps these populations between the vortex pairs, 
the energy behavior of the system will be different when the vortices from a single pair 
are brought close together before and after the braiding. Detecting this energy shift 
reveals the non-Abelian statistics. 

By evaluating the holonomy, we were able to identify a range of the model's 
parameters where the simulation approximates well the exact time evolution of a physical 
system and where the statistics corresponds to the Ising anyons. As expected, larger 
systems exhibit the predicted statistics with higher fidelity. The required magnitude of 
K, however, is larger than anticipated (K 0.1 for parametrization (iii)). In the original 
work [9], the three-body term appears in third order perturbation theory when one 
considers a general Zeeman term (h of) as a perturbation. In our normalization 

the expansion is valid when h 2 <C 1. Since K ~ h 3 , for K = 0.1 one can estimate 
h 2 w 0.2, which clearly does not satisfy the criteria. Therefore, in order to introduce 
the three-body terms into the Hamiltonian perturbatively, such as by adding a small 
magnetic field in the optical lattice proposal [15J, one needs to consider larger systems. 
On the other hand, were the three-body terms engineered [17], our calculation provides 
exact predictions for braiding experiments in such systems. 

5. Conclusions 

In summary, we formulated a method to directly study non-Abelian statistics in 
exactly solvable lattice models whose ground state admits representation as a Slater 
determinant. By applying it to Kitaev's model, we identified finite regions of the 
couplings, where the non-Abelian statistics corresponds to Ising anyons. This confirms 
the previous conjectures for the presence of a non-Abelian topological phase. Finally, 
we proposed a scheme for the implementation and detection of non-Abelian statistics in 
the laboratory. Such an experiment would be an important step towards the physical 
realization of topological quantum computation. It is an interesting topic for future 
research to study whether the holonomy can be used as an order parameter for the 
topological phase when the system is subject to perturbations. 
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